/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TASMatrix
  \brief This class calculates S-matrix exp(i*chi), from given eikonal phase chi.
  Mostly, this class fits the S-matrix with sum of Gaussians, and output the
  fitting result alphaj, i.e., exp(i*chi)=\sum_j{(alphajR+i*alphajI)*exp(-b^2/betaj)},
  where betaj=RL/j.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/08
  \date Last modified: 2020/09/06 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#ifndef _TASMatrix_h_
#define _TASMatrix_h_

#include <complex>
#include <string>
#include "TAException.h"

using cdouble = std::complex<double>;
using std::string;

class TAEikonalPhase;

class TASMatrix{
public:
  // densP and densT are nucleonic density distribution of the projectile and the target
  // density = "" means that the density is delta function, i.e. point-like particle
  // usually for nucleons and particles smaller than that
  // ap, at: proj(targ) nucleon size factor: rho_n(r)=exp(-r^2/alphap^2)
  TASMatrix(int zP, int aP, int zT, int aT, double ek, const string &densP,
    const string &densT, double ap, double at);
  // opt: the optical potential
  TASMatrix(int zP, int aP, int zT, int aT, double ek, const string &opt,
    double ap, double at);
  virtual ~TASMatrix();

  /// calculate and output the S-matrix
  double SMatrixR(double b) const;
  double SMatrixI(double b) const;
  double SMatrix2(double b) const;
  cdouble SMatrix(double b) const;
  double PM() const; ///< projectile mass in MeV
  double TM() const; ///< target mass in MeV
  /// real and imag part of S-matrix, in an array (0, kBMax) with step h=kBMax/fNB
  void SMatrix(double *b, double *smr, double *smi, double *sm2);
  /// fit the S-matrix to an expansion of Gaussians
  /// the results are stored in input arrays alphaj and betaj and members fAlphajR,I
  void GaussianFit();
  TAEikonalPhase *EikonalPhase() const{ return fEikonalPhase; }

  double *GetAlphajR();
  double *GetAlphajI();
  double GetRL() const{ return fRL; }
  void SetRL(double rl){ fRL = rl; }
  int GetNGaus() const{ return fNG; }
  double BMax() const{ return fBMax; }
  /// gaussians used for S-Matrix expansion, \retval array funci[ma] evaluated at x
  static void gaus(double x, double *funci, int ma, const double *p);

protected:
  TAEikonalPhase *fEikonalPhase; ///< contains everything needed to calculate eikonal phase
  int fNB; ///< number of b sampling for eikonal phase
  int fNG; ///< number of Gaussians used in the expansion
  /// beta[j] = fRL/j; fRL is typical of nuclear size, e.g. RL~50fm for 15C+B9
  double fRL;
  double *fB, *fSmr, *fSmi, *fSm2; // Sm(b) = Smr(b) +iSmi(b)
  /// the fitting results
  double *fAlphajR, *fAlphajI;
  /// as the name indicates, unit: fm
  double fBMax;
  bool fFitted; // whether GaussianFit() has been called or not
  bool fSmxed; // if s-matrix has been calculated
};

#endif
